function IRF = BVARir_single(VAR, sigma_draw, beta_draw, impulse, nsteps)
% =======================================================================
% Compute the Choleski IRFs for a VAR model with a given VCV matrix and a
% given matrix of estimated parameters. ONLY THE ONE shock specified in
% impulse is computed. This program can be used to get the IR of a VAR 
% estimated with Bayesian methods (with BVAR_drawpost.m for example). 
% =======================================================================
% IRF = BVARir(VAR, sigma_draw, beta_draw, impulse, nsteps)
% -----------------------------------------------------------------------
% INPUT
%   VAR        : structure, result of VARmodel function
%   sigma_draw : draw of VCV matrix of residuals (obtained with BVARpostdraw or original Choleski)
%   beta_draw  : draw of beta coeff of VAR (obtained with BVARpostdraw or original Choleski)
%   impulse    : vector specifying the position of the variable to shock
%   nsteps    : number of nsteps to compute the IRFs
%
% OUPUT
%   IRF(t,j)   : matrix with the IRFs (t steps, j variable)
% =======================================================================
% Ambrogio Cesa Bianchi, May 2012
% ambrogio.cesabianchi@gmail.com



%% Retrieve parameters and preallocate variables
%===============================================
c_case   = VAR.c_case;
nvars    = VAR.neqs;
nlags    = VAR.nlag;
low_chol = chol(sigma_draw)';


%% Compute the impulse response
%==============================
% check for constant, trend, and/or exogenous variables and remove them
switch c_case
    case 0
        beta_draw_t = beta_draw';
        big_gamma = [beta_draw_t(:,1:nvars*nlags)
            eye(nvars*(nlags-1)) zeros(nvars*(nlags-1),nvars)];

    case 1
        beta_draw_t = beta_draw';
        big_gamma = [beta_draw_t(:,2:nvars*nlags+1)
            eye(nvars*(nlags-1)) zeros(nvars*(nlags-1),nvars)];

    case 2
        beta_draw_t = beta_draw';
        big_gamma = [beta_draw_t(:,3:nvars*nlags+2)
            eye(nvars*(nlags-1)) zeros(nvars*(nlags-1),nvars)];
end

big_gamma_pot = eye(size(big_gamma,1));

% Initialize the impulse response vector
imp_res = zeros(nvars, nsteps);

% First period impulse response (=impulse vector)
imp_res(:,1) = low_chol*impulse;
big_impulse  = [(low_chol*impulse)' zeros(1, nvars*(nlags-1))]';

% Recursive computation of impulse response
for kk = 2:nsteps
    big_gamma_pot = big_gamma_pot * big_gamma;
    big_imp_res   = big_gamma_pot * big_impulse;
    imp_res(:,kk) = big_imp_res(1:nvars);
end
IRF = imp_res';
